# Rescaling coal marginal damages
library(pacman)
p_load(here, fst, data.table,ggplot2,scales)

# NOTE: This script doesnt work any more by itself 
# need tot_md_dt_clean_coal from 05-marginal-damages-oge.R if you want
# to run it 

# Checking results ------------------------------------------------------------
tot_md_dt_clean_coal[,.(
    median = median(value),
    sd = sd(value),
    median_adj = median(value_adj),
    sd_adj = sd(value_adj)),
    keyby = .(fc = fuel_category)
]
# Histogram of damages
ggplot(
    tot_md_dt_clean_coal, 
    aes(
        x = value_adj/1000, # Switch to value to see original
        fill = ifelse(fuel_category %in% c('coal','natural_gas'),fuel_category, 'other')
    )
) + 
geom_histogram( 
    alpha= 0.5, 
    bins = 100,
    position = 'Identity'
  ) +
  scale_x_log10(labels = label_dollar(accuracy = 0.001)) + 
  theme_minimal() + 
  labs(
    x = 'Marginal Damage ($/KWh)',  
    y = 'Number of Power Plants'
  )+   
  scale_fill_brewer(
    name = 'Fuel Category',
    palette = 'Dark2'
  ) +
  facet_wrap(
    ~ifelse(variable == 'md_per_mwh_high', 'Above median production','Below median production'),
    scales = 'free_y',
    nrow = 2, ncol = 1
  ) + 
  theme(legend.position = 'bottom')


ggplot(tot_md_dt_clean_coal[fuel_category == 'coal']) + 
    geom_density(aes(x = value/1000)) + 
    geom_density(aes(x = value_adj/1000), color = 'red') + 
    scale_x_log10() 
